#include <emmintrin.h>

void sgemm( int m, int n, float *A, float *C )
{

  __m128 j1,j2,j3,j4,c1,i_col1, i_col2, i_col3, i_col4, i_col5;

  int i4, i8, i12, i16;
  int top_m = m/20*20;
  int top_m2 = m/4*4;
  for( int k = 0; k < n; k++ )
  {
    
    for( int j = 0; j < top_m2; j+=4 )
    {
      int stored = j+1, stored2=j+2, stored3=j+3;
      j1 = _mm_load1_ps(A+(j+k*m));
      j2 = _mm_load1_ps(A+((stored)+k*m));
      j3 = _mm_load1_ps(A+((stored2)+k*m));
      j4 = _mm_load1_ps(A+((stored3)+k*m));
      for( int i = 0; i < top_m; i+=20 )
      {
	i4=i+4;
	i8=i+8;
	i12=i+12;
	i16=i+16;
        //C[i+j*m] += A[i+k*m] * A[j+k*m];
        c1 = _mm_loadu_ps(C+(i+j*m));
        i_col1 = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j1));
        _mm_storeu_ps(C+(i+j*m),c1);

        c1 = _mm_loadu_ps(C+((i4)+j*m));
        i_col2 = _mm_loadu_ps(A+((i4)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col2,j1));
        _mm_storeu_ps(C+((i4)+j*m),c1);

        c1 = _mm_loadu_ps(C+((i8)+j*m));
        i_col3 = _mm_loadu_ps(A+((i8)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col3,j1));
        _mm_storeu_ps(C+((i8)+j*m),c1);

        c1 = _mm_loadu_ps(C+((i12)+j*m));
        i_col4 = _mm_loadu_ps(A+((i12)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col4,j1));
        _mm_storeu_ps(C+((i12)+j*m),c1);

        c1 = _mm_loadu_ps(C+((i16)+j*m));
        i_col5 = _mm_loadu_ps(A+((i16)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col5,j1));
        _mm_storeu_ps(C+((i16)+j*m),c1);

        //unroll j
	c1 = _mm_loadu_ps(C+(i+(stored)*m));
        //i_col = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j2));
        _mm_storeu_ps(C+(i+(stored)*m),c1);

        c1 = _mm_loadu_ps(C+((i4)+(stored)*m));
        //i_col = _mm_loadu_ps(A+((i4)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col2,j2));
        _mm_storeu_ps(C+((i4)+(stored)*m),c1);

        c1 = _mm_loadu_ps(C+((i8)+(stored)*m));
        //i_col = _mm_loadu_ps(A+((i8)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col3,j2));
        _mm_storeu_ps(C+((i8)+(stored)*m),c1);

        c1 = _mm_loadu_ps(C+((i12)+(stored)*m));
        //i_col = _mm_loadu_ps(A+((i12)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col4,j2));
        _mm_storeu_ps(C+((i12)+(stored)*m),c1);

        c1 = _mm_loadu_ps(C+((i16)+(stored)*m));
        //i_col = _mm_loadu_ps(A+((i16)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col5,j2));
        _mm_storeu_ps(C+((i16)+(stored)*m),c1); 

	c1 = _mm_loadu_ps(C+(i+(stored2)*m));
        //i_col = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j3));
        _mm_storeu_ps(C+(i+(stored2)*m),c1);

        c1 = _mm_loadu_ps(C+((i4)+(stored2)*m));
        //i_col = _mm_loadu_ps(A+((i4)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col2,j3));
        _mm_storeu_ps(C+((i4)+(stored2)*m),c1);

        c1 = _mm_loadu_ps(C+((i8)+(stored2)*m));
        //i_col = _mm_loadu_ps(A+((i8)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col3,j3));
        _mm_storeu_ps(C+((i8)+(stored2)*m),c1);

        c1 = _mm_loadu_ps(C+((i12)+(stored2)*m));
        //i_col = _mm_loadu_ps(A+((i12)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col4,j3));
        _mm_storeu_ps(C+((i12)+(stored2)*m),c1);

        c1 = _mm_loadu_ps(C+((i16)+(stored2)*m));
        //i_col = _mm_loadu_ps(A+((i16)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col5,j3));
        _mm_storeu_ps(C+((i16)+(stored2)*m),c1);

	c1 = _mm_loadu_ps(C+(i+(stored3)*m));
        //i_col = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j4));
        _mm_storeu_ps(C+(i+(stored3)*m),c1);

        c1 = _mm_loadu_ps(C+((i4)+(stored3)*m));
        //i_col = _mm_loadu_ps(A+((i4)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col2,j4));
        _mm_storeu_ps(C+((i4)+(stored3)*m),c1);

        c1 = _mm_loadu_ps(C+((i8)+(stored3)*m));
        //i_col = _mm_loadu_ps(A+((i8)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col3,j4));
        _mm_storeu_ps(C+((i8)+(stored3)*m),c1);

        c1 = _mm_loadu_ps(C+((i12)+(stored3)*m));
        //i_col = _mm_loadu_ps(A+((i12)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col4,j4));
        _mm_storeu_ps(C+((i12)+(stored3)*m),c1);

        c1 = _mm_loadu_ps(C+((i16)+(stored3)*m));
        //i_col = _mm_loadu_ps(A+((i16)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col5,j4));
        _mm_storeu_ps(C+((i16)+(stored3)*m),c1);  

      }
    }
  }

  if(m%4!=0)
  {
  int temp_m = m-1;
  for(int k=0;k<n;k++)
  {
    for(int j=(m-(m%4)); j<m;j++)
    {
    j1 = _mm_load1_ps(A+(j+k*m));
    for( int i = 0; i < top_m; i+=20 )
      {
	i4=i+4;
	i8=i+8;
	i12=i+12;
	i16=i+16;
        //C[i+j*m] += A[i+k*m] * A[j+k*m];
        c1 = _mm_loadu_ps(C+(i+j*m));
        i_col1 = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j1));
        _mm_storeu_ps(C+(i+j*m),c1);

        c1 = _mm_loadu_ps(C+((i4)+j*m));
        i_col1 = _mm_loadu_ps(A+((i4)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j1));
        _mm_storeu_ps(C+((i4)+j*m),c1);

        c1 = _mm_loadu_ps(C+((i8)+j*m));
        i_col1 = _mm_loadu_ps(A+((i8)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j1));
        _mm_storeu_ps(C+((i8)+j*m),c1);

        c1 = _mm_loadu_ps(C+((i12)+j*m));
        i_col1 = _mm_loadu_ps(A+((i12)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j1));
        _mm_storeu_ps(C+((i12)+j*m),c1);

        c1 = _mm_loadu_ps(C+((i16)+j*m));
        i_col1 = _mm_loadu_ps(A+((i16)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j1));
        _mm_storeu_ps(C+((i16)+j*m),c1);

      }
      }
    }
  }

  if(m%20>0)
  {
    int boundary20 = (m-(m%20));
    for( int k = 0; k < n; k++ )
    {
    
      for( int j = 0; j < m; j++ )
      {
        j1 = _mm_load1_ps(A+(j+k*m));
        for( int i = boundary20; i < (m-(m%4)); i+=4 )
        { 
        //C[i+j*m] += A[i+k*m] * A[j+k*m];
        c1 = _mm_loadu_ps(C+(i+j*m));
        i_col1 = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j1));
        _mm_storeu_ps(C+(i+j*m),c1);
        }
      }
    }
  }

  if(m%4>0)
  {
    int boundary4 = (m-(m%4));
    for(int k=0; k<n;k++)
      for(int j=0; j<m;j++)
	for(int i=boundary4; i<m;i++)
	  C[i+j*m] += A[i+k*m] * A[j+k*m];
  }
}

